function [A, dAdpar] = T1IRMagnitudeAbs( par, echoes)

% computes the magnitude M with given parameters and inversion times
% echoes = [n x 1] inversion times
% par = the 3 element parametervector [A, B, R1]'
%
% The datamodel is given by
% M = A - B*exp(-echoes*R1).
%
% Created by Henk Smit, EMC, 01-2011 based on the work by Dirk Poot, University of Antwerp, 13-8-2007

if size(par,1) ~= 3 %HENK
    error('incorrect parameter vector');
end;
numtr = size(par,2);
numgr = size(echoes,1);

A = zeros(numgr,numtr);
A_ex = A;
A_part = A;
Aabs = A;

for k=1:numtr
 A_ex(:,k) = exp(-echoes(:,1)*par(3,k)); %HENK matrix with all results 
 A_part(:,k) = par(2,k) * A_ex(:,k);
 Aabs(:,k) = par(1,k) - par(2,k) * A_ex(:,k);   
 A(:,k) = abs(par(1,k) - par(2,k) * A_ex(:,k));   
end;

if nargout>1
    dAdpar=([A(:,1) A(:,1) A(:,1)]./[Aabs(:,1) Aabs(:,1) Aabs(:,1)]).*[ones(size(A_ex,1),size(A_ex,2)) -A_ex(:) repmat(echoes,numtr,1).*A_part(:)]; %HENK 
end;
